2.1.1 Overview

Structures such as a street or a building can be represented by a geometric objects. For instance, a street can be represented by a straight line, whereas a building by a polygon. It makes sense to store the data of these type of structures as their geometric representation. For example, the building we are in and the adjacent street:

# This code in incomprehensible to you FOR NOW
# Just an example of geometries
r <- rast("data/idiv-building.tif")
plotRGB(r)
l <- vect(
  rbind(
    c(12.39585, 51.31866),
    c(12.39846, 51.31722)
  ),
  crs = "EPSG:4326"
)
lines(l, col = "green", lw = 5)
idiv <- vect(
  rbind(
    c(12.39573, 51.31794),
    c(12.39645, 51.31765),
    c(12.39664, 51.31812),
    c(12.39602, 51.31849),
    c(12.39573, 51.31794)
  ),
  crs = "EPSG:4326"
) |>
  as.lines() |>
  as.polygons()

## Warning: [as.polygons] dropped attributes

idiv$name <- "iDiv building"
idiv$perimeter <- perim(idiv)
idiv$area <- expanse(idiv)
plot(idiv, add = TRUE, col = adjustcolor("blue", alpha.f = .5))

Geometries are saved in vector or shape files. In addition to the geometric representation of the structure, shapefiles also contain metadata for the geometric representation, such as the its extent and coordinate reference system (we will talk about this later), and for the structure, such as their name, length, etc. For instance, the geometry of the iDiv building above is:

idiv

##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 1, 3  (geometries, attributes)
##  extent      : 12.39573, 12.39664, 51.31765, 51.31849  (xmin, xmax, ymin, ymax)
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       :          name perimeter  area
##  type        :         <chr>     <num> <num>
##  values      : iDiv building     237.7  3265

Note that this structure has one geometry and three attributes/variables (dimensions : 1, 3 (geometries, attributes)): name, perimeter (\(m\)), and area (\(m^2\)).

2.1.2 Geometry Types

In terra, there are three main geometry types:

1. Points: they are defined by a vector with two values for coordinates (e.g., longitude and latitude)
2. Lines: they are a series of points connected pairwise by straight lines
3. Polygons: they are a series of lines inscribing a closed area

In terra, geometries are created with the function vect(). To create a point geometry, you need to pass a matrix with the coordinates:

xy <- cbind(12.39585, 51.31866)  # cbind force it to be a matrix
vect(xy)

##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 1, 0  (geometries, attributes)
##  extent      : 12.39585, 12.39585, 51.31866, 51.31866  (xmin, xmax, ymin, ymax)
##  coord. ref. :

To create a geometry of multiple points, simply pass a matrix that has multiple rows:

xy <- rbind(
  cbind(12.39585, 51.31866),
  cbind(12.39584, 51.31865)
)
poi <- vect(xy)

In order to create lines, first create a multi-point geometry, then convert it using as.lines():

xy <- rbind(
  cbind(12.39585, 51.31866),
  cbind(12.39584, 51.31865)
)
lin <- vect(xy) |> as.lines()

In order to create polygons, first create a mutli-point geometry, then covert it to lines, and finally to polygons using as.polygons(). The last point should be the same as the first to close the geometry.

xy <- rbind(
  cbind(12.39585, 51.31866),
  cbind(12.39584, 51.31865),
  cbind(12.39584, 51.31866),
  cbind(12.39585, 51.31866)
)
pol <- vect(xy) |> as.lines() |> as.polygons()

## Warning: [as.polygons] dropped attributes

plot(pol, col = "gold")
lines(lin, col = "red", lw = 5)
points(poi, col = "blue", cex = 1)

2.1.3 Reading & Writing Geometries

vect() is also used to read geometries from a file. For instance, to read the shapefile with the country boundary of Germany:

ger <- vect("data/germany.shp")
ger

##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 1, 5  (geometries, attributes)
##  extent      : 5.866755, 15.04179, 47.27012, 55.05866  (xmin, xmax, ymin, ymax)
##  source      : germany.shp
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       :            shapeName shapeISO              shapeID shapeGroup
##  type        :                <chr>    <chr>                <chr>      <chr>
##  values      : the Federal Republi~      DEU 19620994B6459175825~        DEU
##  shapeType
##      <chr>
##       ADM0

To write geometries to disk, use writeVector().

writeVector(ger, "data/germany.shp", overwrite = TRUE)

2.2.1 Overview

Rasters are gridded areas where each grid pixel assumes a value. When the data we are interested in are values on a gridded area, rasters are a better options than to using geometries. For example, if we want to represent landcover type over a large area(forest, agriculture, build-up, etc.), it is easier to grid the area into pixels and save the landcover type of each pixel rather than to create a geometry for each different structure.

r <- rast("data/landcover-2015.tif")
plot(r)  # from EOC of DLR

Rasters are saved in raster files. In addition to the values of each grid, rasterfiles also contain metadata for the grid, such as the its extent and coordinate reference system (we will talk about this later). For instance, some of the metadata of the landcover raster:

r

## class       : SpatRaster 
## dimensions  : 2042, 2032, 1  (nrow, ncol, nlyr)
## resolution  : 10, 10  (x, y)
## extent      : 4477870, 4498190, 3126660, 3147080  (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
## source      : landcover-2015.tif 
## color table : 1 
## categories  : category 
## name        :        category 
## min value   : Artificial Land 
## max value   :           Water

Note the metadata dimensions (grid size), resolution (spatial resolution, in \(m\)), and extent (spatial extent).

2.2.2 Creating Rasters from Scratch

Simply put, rasters are array associated spatial metadata. In terra, use rast() to create a raster from a matrix.

m <- matrix(rnorm(100), nrow = 10, ncol = 10)  # 10x10 matrix
r <- rast(m)
plot(r)

Usually, you will not create rasters from scratch. However, if you are familiar with R arrays, you will notice a nice properties that is inherited by rasters: rasters can have multiple layers, with each layers being a matrix.

r <- rast(array(rnorm(200), dim = c(10, 10, 2)))
plot(r)

When an area has multiple rasters, it is possible to stack them to create a single raster object. The rasters in the stack need to have all the same extent and coordinate reference system. Use c() to stack rasters using terra.

r1 <- rast(matrix(rnorm(100), nrow = 10, ncol = 10))
r2 <- rast(matrix(rnorm(100), nrow = 10, ncol = 10))
r <- c(r1, r2)  # a stack
plot(r)

2.2.3 Reading & Writing Rasters

To read rasters into memory, use rast().

r <- rast("data/landcover-2015.tif")
r

## class       : SpatRaster 
## dimensions  : 2042, 2032, 1  (nrow, ncol, nlyr)
## resolution  : 10, 10  (x, y)
## extent      : 4477870, 4498190, 3126660, 3147080  (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
## source      : landcover-2015.tif 
## color table : 1 
## categories  : category 
## name        :        category 
## min value   : Artificial Land 
## max value   :           Water

To write rasters to disk, use writeRaster().

writeRaster(r, "data/landcover-2015.tif", overwrite = TRUE, datatype = "INT1U")

Here, the argument datatype = "INTU1" (specifying the values in the grid are 1-byte unsigned integers) is needed because the rasters contains classes (levels). In most cases, this argument is best left unspecified, as terra picks the optimal value by default.

2.3.1 Spatial extent

The spatial extent is the quadrilateral that inscribe the area of the spatial data. The extent is usually represented by the coordinates of the four vertices of the quadrilateral, i.e. xmin, xmax, ymin, and ymax. To get the spatial extent of an object, use ext().

r <- rast("data/landcover-2015.tif")
ext(r)

## SpatExtent : 4477870, 4498190, 3126660, 3147080 (xmin, xmax, ymin, ymax)

2.3.2 Coordinate reference system

GIS try to represent the surface of the Earth, a 3D spheroid, onto a plane. The coordinate reference system (CRS), also known as spatial reference system (SRS), defines how this projection of a 3D object to a 2D one is achieved. It is not possible to achieve this projection accurately and some distortions will always be present. In particular, at least one of distance, angular conformity, and area will be distorted. Projections can be grouped into types, depending on which property of the Earth surface they do not distort:

• Conformal projections: they correctly represent the angles between points and, thus, shapes (e.g. ESRI:54004¹).
• Equidistant: they correctly represent distances (e.g. ESRI:54002).
• Equal-area: they correctly represent areas (e.g. ESRI:54034).

An overview of ESRI and EPSG projections can be found at https://spatialreference.org/. Wikipedia also has a nice list with the property of each projection: https://en.wikipedia.org/wiki/List_of_map_projections.

To know the CRS of an object, use crs(). By default, crs() displays the CRS in Well Known Text (WKT) format.

crs(idiv, parse = TRUE)

##  [1] "GEOGCRS[\"WGS 84\","                                   
##  [2] "    DATUM[\"World Geodetic System 1984\","             
##  [3] "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,"   
##  [4] "            LENGTHUNIT[\"metre\",1]]],"                
##  [5] "    PRIMEM[\"Greenwich\",0,"                           
##  [6] "        ANGLEUNIT[\"degree\",0.0174532925199433]],"    
##  [7] "    CS[ellipsoidal,2],"                                
##  [8] "        AXIS[\"geodetic latitude (Lat)\",north,"       
##  [9] "            ORDER[1],"                                 
## [10] "            ANGLEUNIT[\"degree\",0.0174532925199433]],"
## [11] "        AXIS[\"geodetic longitude (Lon)\",east,"       
## [12] "            ORDER[2],"                                 
## [13] "            ANGLEUNIT[\"degree\",0.0174532925199433]],"
## [14] "    USAGE["                                            
## [15] "        SCOPE[\"Horizontal component of 3D system.\"],"
## [16] "        AREA[\"World.\"],"                             
## [17] "        BBOX[-90,-180,90,180]],"                       
## [18] "    ID[\"EPSG\",4326]]"

Notice, among the others, the attributes AREA (the area of usage for the CRS) and ID (in this case, the EPSG code).

WKT is not very nice for humans; use the extra argument proj = TRUE to see the PROJ4 format of the CRS.

crs(idiv, proj = TRUE)

## [1] "+proj=longlat +datum=WGS84 +no_defs"

2.3.3 Resolution

Resolution applies only to rasters, as geometries are geometric representation of structure and can be scaled at any level. To get the resolution of a raster, use res().

r <- rast("data/landcover-2015.tif")
res(r)

## [1] 10 10

The unit of the output of res() is the same as the unit of the CRS, in this case meters:

crs(r, proj = TRUE)

## [1] "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"

5.1.1 Local Summary

Raster stacks can be summarized using common function, such as mean, sd, sum, etc. The output is a single raster layer with the summary. Mean, sum, minimum, and maximum, can be calculated simply using mean(), sum(), etc. For other function, use app().

ndvi <- rast("data/ndvi-2024.tif")
avg_ndvi <- mean(ndvi, na.rm = TRUE)  # equivalent to app(ndvi, "mean")
std_ndvi <- app(ndvi, "sd", na.rm = TRUE)
month_max_ndvi <- app(ndvi, "which.max", na.rm = TRUE)  # month with the maximum NDVI

plot(avg_ndvi, col = colorRampPalette(c("dodgerblue3", "gold", "darkgreen"))(100))

plot(std_ndvi, col = hcl.colors(100, "Reds", rev = TRUE))

This is called a local summary, because it is calculated across pixels with the same coordinates.

5.1.2 Global Summary

To get global summaries, i.e. calculated across all pixels of a layer, use global()

global(ndvi, "mean", na.rm = TRUE)

5.1.3 Zonal Summary

Zonal statistics summarize the values of a raster using categories from another raster of geometry. Use zonal() to obtain zonal statistics.

landcover <- rast("data/landcover-2015.tif")
avg_ndvi <- project(avg_ndvi, landcover)  # this will also resample avg_ndvi (see the chapter Raster Operations)
zonal(avg_ndvi, r, fun = "mean", na.rm = TRUE)

zonal(avg_ndvi, r, fun = "sd", na.rm = TRUE)  # spatial sd

5.2.1 Aggregate

Aggregating creates a raster that is at coarser resolution compared to the original raster. Use aggregate() to aggregate rasters. The arguments fact, number of cells in each direction to be aggregated, and fun, the function used for aggregatation, need to be specified.

r <- rast(matrix(rnorm(16^2), 16, 16))
plot(r)

r_aggr <- aggregate(r, 4, "max")  # maximum values in the 4x4 grid
plot(r_aggr)

5.2.2 Disaggregate

Disaggregating creates a raster that is at finer resolution compared to the original raster. Use disagg() to disaggregate rasters. The arguments fact, number of cells in each direction to create for each cell of the original raster, and method, the method used for disaggregation, need to be specified. Available methods are near, for nearest neighbor, or bilinear, for bilinear interpolation.

r_disagg <- disagg(r, 4, method = "bilinear")
plot(r_disagg)

Aggregating and disaggregating a raster does not give, in general, the same raster back. Some of the original information is lost during aggregation.

plot(c(r, r |> aggregate(4, "mean") |> disagg(4, "bilinear")))